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Abstract — We have developed an approximate signal recovery 
algorithm with low computational cost for compressed sensing on 
the basis of randomly constructed sparse measurement matrices. 
The law of large numbers and the central limit theorem suggest 
that the developed algorithm saturates the Donoho-Tanner weak 
threshold for the perfect recovery when the matrix becomes as 
dense as the signal size N and the number of measurements M 
tends to infinity keep a = A//7V ~ 0(1), which is supported 
by extensive numerical experiments. Even when the numbers of 
non-zero entries per column/row in the measurement matrices 
are limited to 0(1), numerical experiments indicate that the 
algorithm can still typically recover the original signal perfectly 
with an 0(Y) computational cost per update as well if the density 
p of non-zero entries of the signal is lower than a certain critical 
value /3th(a) as A^, M — 5> oo. 

I. Introduction 

Compressed (compressive) sensing (CS) is a framework that 
enables the recovery of a sparse signal from few of its mea- 
surements by exploiting the sparsity as the prior knowledge 
of the original signal. Famous applications of CS include 
computational photography and seismic data processing, in 
which the original signal and its measurements are linearly 
related, in most cases, by the Fourier and/or wavelet trans- 
formations. In general, these relationships are expressed by 
dense matrices. Accordingly, much effort has been made in 
analyzing the performance IT], (IS, iB] and developing 
practical algorithms for the signal recovery ||5], (0, (|7] for 
the density matrix based CS. 

Data compression, data stream processing, and group testing 
are also included in examples of CS but may be classified 
into another category. Unlike the above applications, these 
examples naturally allow situations where each measurement 
is related to only a few entries of the original signal. In such 
cases, it is more suitable to model the signal-measurements 
relationship by not density but sparse matrices [8|. 

The Zi-norm minimization under the constraint of given 
measurements is widely adopted for the signal recovery of 
CS and the interior point method is a standard solution for 
this task. In general, the computational cost per update of 
this scheme enlarges as cubic of the problem size, which 
is regarded as feasible in computational complexity theory. 
However, in practice, the cost may exceed the allowable limit 
depending on the problem size and/or situations. In such cases. 



the sparsity of the matrices can be advantageous for reducing 
the computational cost because necessary computations for 
relating the signal and measurements increase just linearly 
with respect to the signal size if the numbers of non-zero 
entries per column/row in the matrix are 0(1) |l8l, |j9l- 

From this perspective, we have developed an approximate 
algorithm for the signal recovery of the sparse matrix based 
CS. In a similar setting, an earlier study llTOl developed an 
algorithm on the basis of belief propagation (BP) fTTl, fT2l 
in conjunction with using a mixture of two finite variance 
Gaussians as a sparse prior The computational cost of the 
algorithm increases only linearly with respect to the signal 
size. However, the following two issues may be problematic: 

• Functions of continuous variables must be updated to deal 
with continuous signals. This means that the prefactor 
of the computational cost is rather large even if it is 
proportional to the signal size. 

• The finiteness of the variances of the two Gaussians in the 
prior does not allow the perfect recovery of the original 
signal even if no noise is added to the measurements. 

We will show that these possible drawbacks can be resolved 
by using BP to the constrained Zi-norm minimization in 
conjunction with further quadratic approximation. We will 
also show that our algorithm asymptotically saturates the 
Donoho-Tanner weak threshold for the perfect recovery [1T31 
when the matrix becomes dense. This indicates that, in the 
dense limit, the algorithm developed here is as capable as the 
approximate message passing (AMP) proposed by Donoho et 
al. Q although control of a soft-thresholding parameter is not 
necessary unlike AMP. 

The remainder of this article is organized as follows. The 
next section introduces the problem on which we will focus. 
Section III, which is the main part of this article, describes the 
suitable signal recovery algorithm for the sparse matrix based 
CS that we developed. Section IV examines the properties of 
the algorithm by extensive numerical experiments. The final 
section summaries the paper. 

II. The problem setting 

In the general scenario, we suppose that each entry e M 
of an iV-dimensional original signal = {xD E is 
independently generated from an identical distribution P{x) = 



(1 — p)S{x) + pP{x). Here, P{x) is a certain distribution 
that has a finite variance but not a finite mass at the origin, 
and p represents the density of non-zero entries of the signal. 
Multiplying M{< N)xN matrix F = (F^,) e M*^^^ yields 
the linear measurements y = Fx'^ of dimensionality M. To 
simply characterize the sparsity of the measurement matrix, 
we focus on the case of regular ensembles in which positions 
of non-zero entries in F are determined randomly under the 
constraint that the numbers of non-zero entries per column 
and row are fixed to finite values j and k, respectively while 
N,M -> oo with keeping a — M/N ^ 0{l). However, 
extension to irregular ensembles for which values of j and 
k are distributed in a matrix is straightforward. After the 
positions of the non-zero entries are fixed, each value is 
provided by independently sampling a random number from 
an identical distribution of a finite variance. 

To recover the original signal given y and F, we follow 
the constrained 1 1 -norm minimization scheme 

minimize < ^ \xi\ > subject to Fx = y. (1) 

This can be expressed as a problem of linear programming. In 
general, necessary computational cost for solving this scales 
as 0{N^) per update by using the interior point method as 
a certain matrix inversion is required for each iteration [.14l . 
The purpose of our study is to reduce this cost by developing 
an efficient approximate algorithm utilizing the framework of 
BP. 

III. Algorithm development 

A. Belief propagation 

As a basis of our study, we first convert eq. dl} to an 
unconstrained optimization problem 

maxmin |A'^(Fa; - y) + ^ |.Ti|| (2) 

by introducing the Lagrange multiplier A = (A^) G R^^, 
where maxx and mmy stand for maximization and mini- 
mization with respect to X and Y, respectively. T denotes the 
matrix transport. 

The objective function of eq. (|2]l includes non-trivial in- 
teraction terms X^Fx, which is why computation cost is 
so considerable. Our key idea for reducing the cost is to 
approximate eq. (I2]i by a bunch of single-body optimization 
problems that can be handled with a lower computational cost, 
which is similar to the spirit of mean field approximations of 
statistical mechanics 1151 . 

For this, we introduce auxiliary functions 4)i^^{xi) and 
tpfj.^i{^ti) to two types of optimization variables Xi and A^, 
respectively. (t)i^^{xi) physically means the single-body ob- 
jective function of Xi for the "/i-cavity system" that is defined 
by removing A^i from the original system, and similarly for 
fpn^ii^fj.) lT6l . IITtI . The BP framework indicates that these 
functions can be determined by the following considerations: 



1) i-cavity — > /i-cavity 

Recover the original system by inserting Xi to the i- 
cavity system. After that, choose an index p e A4{i), 
where denotes the set of indices of the Lagrange 

multipliers that are directly connected to Xi, and remove 
A^. This yields the /i-cavity system. Optimizing the ap- 
proximate objective function of the /i-cavity system with 
respect to \ueM{i)\ti offers the single body objective 
function of Xi in the /i-cavity system. Here, A\a denotes 
the set that is defined by removing an element a from a 
set A. This provides 

4>l^pi{,Xi) ^\xi\ 

+ ^ TCLa.yi{Fi,i\i,Xi + -ipij^i{Xu)} -0) 

2) /i-cavity i-cavity 

Recover the original system by inserting A^ to the p- 
cavity system. After that, choose an index i E I{p), 
where I{p) denotes the set of indices of the signal 
variables that are directly connected to A^, and remove 
Xi. This yields the i-cavity system. Optimizing the 
approximate objective function of the i-cavity system 
with respect to xi^x{fi)\i offers the single body objective 
function of p in the i-cavity system. This provides 

V^p^jCA^) = -2/^A^ 

+ ^ mm {Fi_aXf_,xi + (t>i^f_,{xi)} . (4) 
iex(ij.)\i ' 

Let us denote the single body objective functions for 
approximating eq. (|2]l as 4'i{xi) and ip^{X^). These can be 
constructed from the auxiliary functions as follows; 

3) Construction of (j)i{xi) and ipfi{Xfi) 

Recover the original system by inserting A^ to the p- 
cavity system and optimize the approximate objective 
function with respect to A^£7vi(j). This yields 

<l)i{Xi)^\xt\+ y2 ma.x{Ffj.iXfj.Xi+ij^^t{X^)} . (5) 
Similarly, 

tpf,{Xf,) = -yf,Xf,+ '^mm{Ff,iXf,Xi + (l)i^f,{xi)} . (6) 

Under suitable conditions, iterating eqs. (|3} and (|4) for all 
connected pairs of i and p is expected to yield a convergent 
solution of (t)i^^{xi) and ■(/;^^i(A^). Inserting the solution 
into (|5]l and optimizing (t>i{xi) offer the recovered signal Xi ~ 

aigmm^. {(j)i{xi)}. 

B. Quadratic approximation 

As long as j, fc ^ 0(1), the necessary computational cost 
for performing the above procedure grows linearly with N 
since the optimization required at each step is concerned with 
only a small number of variables. However, the prefactor is 
considerably large since one has to update functions at each 
step, which may reduce practical applicability. This difficulty 



is shared with another BP-based algorithm developed by Baron 
et al. llOl. 

To reduce necessary computation cost, we limit (f)i^^{xi) 
and ip^^i{X^) to the form of (piecewise) quadratic functions 
as 



of Ai^f^, Bi^^, C/^i^i and D^^i. In addition, substituting eq. 
^ into eq. Q offers an expression 

1 



where 



(7) 



and derive update rules for Ai^^, Bi^f^, Cf^^i and D^^i. 
For this, we first substitute eq. (|8]l into eq. Q, which yields 
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(9) 



(10) 



Equation ( fTSI l provides the recovered signal as 

X, = !{B,-Ai). 



(16) 



(17) 



(18) 



(19) 



exactly. Next, we insert eq. (|7|i into eq. (|4|i. In this, the 
expression 

uAn{F^i\^xi + (j)^,^i{xi)} = 

2A;:r^ ' Bi^^-\^,t^,i>i, 

0, 



All this constitutes the main achievement of this article. Actual 
time required for running this algorithm is much less than that 
for the direct product from BP of eqs. (|3), (E) and (|5) because 
one only has to deal with four types of variables defined 
for each pair of signal variables and the Lagrange multipliers 
without handling their functions, although the growth rate of 
the computational cost with respect to N is the same between 
the two algorithms. 



Bi^fj, A^i^^;|<l, (11) Dense matrix limit 



B, 



KFal<-l 



should be paid attention to. For a fixed /i, eq. (fTTl i represents a 
piecewise quadratic function that switches its functional form 
at two points = [Bi^^ ± 1)/F^i, which vary among I E 
This means that directly using eq. (fTTT l in assessing 
eq. & yields a piecewise quadratic function that switches its 
functional form at 2{j — 1) different points, which makes it 
impossible to obtain the quadratic form of eq. dHJ any more. 
To practically resolve this problem, we approximate the right 
hand side of eq. ( fTTT ) by 

0, \Bi^^ < 1, (12) 



To theoretically examine the property of the algorithm 
developed above, let us suppose an extreme situation where 
j — > M, k ^ N and matrix entries Ff^i independently 
follow an identical distribution of zero mean and variance 
N^^. A consideration similar to the following has been 
provided in research on CDMA multiuser detection of wireless 
communication before [18|. In the supposed situation, the law 



of large numbers suggests that Ai 



A,. 



A EE a/C 



and C^^, ^C = N-^ Ef=iid/dBi)f{Bi; A) typically hold 
for i = 1,2, ... ,N and jj, ~ 1,2, . . . , M, where we used an 
expression of the Taylor expansion 

df{Bi;A)F^i 



f{Bi^^;A)-f{Bi;A) 



B, 



< -1 



(20) 



which is justified if ^ 1 holds. As eq. (fT2] i 

represents not piecewise but totally quadratic functions of A^, 
using this approximation in assessing eq. (|4]i yields the form 
of eq. (|8]l, which leads to 

C,^.= F^^M^t^il^, (13) 



dB, 



D^^,= ^ F^if{Bi^^;Ai^,,). (14) 
Here, f{B;A) represents a soft-thresholding function 



f{B;A)^[B-^]Qm-l)IA, 



(15) 



in which <d{u) ~ 1 for u > and vanishes, otherwise. 

Under appropriate conditions, iterating eqs. (|9) and (ITOi — !• 
eqs. ( fTsT i and ( fT4] i is expected to yield a convergent solution 



- 0(iV-i/2) ^ {N ~^oo). 

Further, inserting this and the expression of j/^ — Yl!i=i Fvi^l 
into eq. ( fTOl i, we have 

B,^, ^^xU^Y. E - /(^'-- ^)) • ^21) 

Concerning this, the central limit theorem indicates that the 
second term of eq. (I2TI 1 converges to obey a Gaussian distri- 
bution of zero mean and variance C~'^N^'^ '^i^^t^ '^i^ii^'i ~ 
f{Bi^.;A))^ ^ («/C2)iV-iE^iK - f{Bf,A))^ and is 
independent among different i's for given ^ as N tends to 
infinity since F^i's are independent of one another. 

These arguments mean that as N, M = aN — 
00, macroscopic variables with respect to the signal es- 
timate of eq. (fT9] l, m = ^^^J2iLi^i^i Q ~ 
Si^i(^j)^ ^re determined from the following equations: 
m = {jDzx'^f{B;A))^,, Q = (jDz{f{B;A)f 



and C = {J Dz{df {B;A) /dB))^„, where B ^ 
{a/C)x° + {y^a{Q~ 2m + Qo)/C)z, A = a/C and Dz = 
dz exp(— z^/2)/V27r denotes the Gaussian measure. (■ • '/^o 
represents the average operation with respect to the original 
signal and Qo ~ ((a;'^)^)^^. It may be noteworthy that 
these equations for determining the macroscopic variables are 
equivalent to those obtained by the replica method of statistical 
mechanics for the /i-norm based signal recovery scheme 14|. 
As the replica method reproduces the result identical to that 
obtained by mathematically rigorous analyses lfT3l . lfT9l . this 
suggests that the current algorithm can saturate a theoretical 
limit of the ^i-norm based scheme for the perfect recovery, 
which is often termed the Donoho-Tanner weak threshold lfT3l , 
in the limit of dense matrices. 

The techniques based on the law of large numbers and the 
Taylor expansion, which are used above, are also useful for 
reducing the necessary computational cost in the dense matrix 
case. Inserting eq. into Z)^ = J^iLi FtJ.if{Bi^ti; A) 

yields 

= + ^ (^^^^^^) (22) 

i— 1 \ i—1 ^ / 

where z^ = y^~ D^, and Yln=i f {Bt^t,\ A) / dB,^^ ~ 
[l/N) df{B,- A)/dB, was used in the right hand side 
on the basis of the law of large numbers. On the other hand, 
plugging D^^j = Df,- Ffj,J{Bi^i,;A) ~ I?^ - Ff^.x, into 
eq. ( fTsT i and using an effect of the law of large numbers 
E*iiF2,;^aleadto 

-Bj = — ^ Ff,iZ^, + -^^i- (23) 

Iterating eqs. (|22] | and ( |23] l in conjunction with C = 
(l/N) J2f=i df{B,;A)/dB, and A = a/C offers the signal 
estimate Xi = f{Bi]A) with an 0{MN) computational cost 
per update, which is considerably smaller than that of the 
original expression, 0{MN{M + N)), when TV, M > 1. 

In the context of the signal recovery of CS, a class of 
algorithms similar to eqs. (l22t and (l2?t . termed AMP, has 
already been proposed by Donoho et al. Q. Performance 
of AMP generally depends on how a certain parameter for 
soft-thresholding, which corresponds to in the current 

algorithm, is controlled externally. A characteristic feature of 
the current algorithm is that this parameter is tuned adaptively 
for given F and y in the course of updates. This property may 
be preferred in practical usage because one does not have to 
care about effects of sample fluctuations in such an adaptive 
control. 

IV. Experimental Validation 

To examine the abilities and limitations of the above 
scheme, we assessed the capability of the perfect signal 
recovery for cases (A) {j,k) — (10,20) and (B) the dense 
matrix case by extensive numerical experiments. As they 
can be compared with accurate and mathematically rigorous 
assessments, we will first show the results for (B). 



In the (B) experiment, the probabiUty that the original signal 
is perfectly recovered up to 10000 iterations of eq. (l22i — > 
C = {l/N)j:f^^dfiB,;A)/dB, ^ A = a/C ^ eq. m 

eq. fT% was evaluated through 10000 trials for each pair of 
various signal density p and signal length N. We focused on 
the case of a = M/N — 1/2. For each trial, we generated an 
M(= A^/2) xiV dense random matrix F, entries of which were 
sampled independently from an Gaussian distribution of zero 
mean and variance 1/iV. Each entry of the original signal 
x'-* was sampled from P{x) = (1 — p)S{x) + pexp (— a;^/2) 
independently. We judged that x'^ was perfectly recovered 
if N-^J2Ziixt - xff < 10-^ is satisfied. The results 
for N =500, 1000, and 2000 are plotted in Fig. [T] which 
shows that three curves of the probability of the perfect 
recovery for the different system sizes intersect one another 
at a point very close to the Donoho-Tanner weak threshold 
for a = 1/2, Pc(l/2) = 0.1928.... This suggests that the 
developed algorithm can saturate the theoretical limit of the 
perfect recovery in the dense matrix case with a computational 
cost of 0(7V^) per update, which is lower than that required 
for the generic interior point method. 

For the case of (A), we also assessed the probability of the 
perfect recovery by using iterations of eqs. (|9) and (ITffi 
eqs. ( fTsT i and (fT4l l for sparse matrices which were randomly 
constructed under the constraint of (j, k) — (10, 20). Values of 
non-zero entries were independently sampled from the Gaus- 
sian of zero mean and unit variance. As convergence is rather 
faster than in the dense matrix case, we set the maximum 
number of iterations to 1000. The ways of generating a;", the 
condition for judging the perfect recovery, and the number of 
trials for each parameter setting were the same as the above. 
Figure |2] (a) plots the results for N =3200, 6400, 12800, and 
25600 and shows that four curves for the different system sizes 
non-trivially intersect one another at p ^ 0.1652. The accurate 
estimate for the theoretical limit of the perfect recovery 
has not been clarified for the sparse matrices. Alternatively, 
we performed another experiment by shuffling connectivities 
among variables at each iteration, which corresponds to the 
density evolution analysis of BP l.20il . The results are presented 
in Fig.[2](b). A non-trivial cross of four curves is also observed 
at p ^ 0.1643, where the difference in the third digit from that 
in Fig. |2] (a) may be attributed to effects of finiteness of the 
system sizes and the maximum number of iterations. These 
suggest that for the sparse matrices the algorithm developed 
here can typically recover the original signal perfectly with an 
0{Nk + Mj) cost of computations per update if the density of 
non-zero entries in the signal is below a certain finite critical 
value pth(a) as N,M ~^ oo keeping a — M/N ~ 0(1). 

V. Summary 

In summary, we have developed an approximate signal 
recovery algorithm with low computational cost for sparse 
matrix based compressed sensing. The developed algorithm 
saturates the Donoho-Tanner weak threshold for the perfect 
recovery when the measurement matrix becomes dense. For 
sparse matrices, the algorithm is still capable of typically 
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Fig. 1. Probability of perfectly recovering the original signal for dense 
matrices. 
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Fig. 2. (a): Probability of perfectly recovering the original signal evaluated 
for randomly constructed sparse matrices of (j, fc) = (10, 20). (b): That 
assessed by the "density evolution" scheme (see the main text for details). 



recovering the original signal perfectly if the density of the 
non-zero entries in the original signal is lower than a certain 
finite critical value as the signal length and the number of 
measurements tends to infinity while keeping the ratio between 
them finite. 

Accurately identifying the theoretical limit of the perfect 
signal recovery for the sparse matrix based compress sensing 
and designing sparse matrices fJTl for improving the perfor- 



mance of the current algorithm are challenging problems for 
future studies. 
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